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Abstract. We study the non-equilibrium time evolution of the classical XY spin model in two dimensions. 
The two-time autocorrelation and linear response functions are considered for systems initially prepared 
in a high temperature state and in a completely ordered state. After a quench into the critical phase, we 
determine, via Monte Carlo simulations, the time-evolution of these quantities and extract the temperature 
dependence of the slope of the parametric plot susceptibility/correlation in the asymptotic regime. This 
slope is usually identified with the infinite fluctuation-dissipation ratio which measures the violation to the 
equilibrium fluctuation-dissipation theorem. However, a direct measure of this ratio leads to a vanishing 
value. 

PACS. 75.40.Gb Dynamic properties 

05.70.Ln Non-equilibrium and irreversible thermodynamics 



1 Introduction 

Non-equilibrium properties of classical spin systems have 
received a lot of interest these last years, especially 
in the context of aging 1 , 2 and from the point of 
view of fluctuation-dissipation theorem (FDT) and its 
extensions. |H] The main feature of non-equilibrium dy- 
namics is the breakdown of time-translation invariance, 
which is a characteristic that has been used recently, 
together with space-symmetries in order to build a 
space-time conformal like theory for some scale invariant 
systems. 0] This theory has given some predictions that 
have been already tested on some systems, like the Ising 
with Glauber dynamics or spherical model. [5] A system 
relaxing towards its equilibrium state shows, in the ag- 
ing regime, a dependence in the two-time functions on 
both the observation time t and the so called waiting time 
t w < t. In this context, an extension of the fluctuation- 
dissipation theorem (FDT) was proposed. At equilib- 
rium, the FDT relates the correlation function to its con- 
jugate linear response function via 



It measures the ratio between the actual response and the 
expected response if the FDT was valid. Recently, a lot of 
interest was put in the asymptotic value of the FDT ratio, 
defined by 



X- 



lim lim X(t, t w ) 

t w — >oo t — >oo 



(3) 



In particular, Godreche and LuckJTj proposed that this 
quantity should be universal for a critical quench. Evi- 
dences to support this universality were obtained on ex- 
actly solvable spin systems quenched from infinite tem- 
perature, numerically on 2d and 3d Ising model with 
Glauber dynamics and also checked from field-theoretic 
two-loop expansion of the 0(n) model 9 and the 2d voter 
model. [lOj This universality was tested recently for a wide 
class of initial states in Id Glauber Ising model. [TT] 

However, going on more complex systems the linear 
response function itself is not accessible by numerical, nor 
experimental analyses and one is forced to look on inte- 
grated response functions, that is measuring susceptibili- 
ties: 



R(t-t w ) = P^-C{t-t w ) , 
ot w 



(1) 



X (t,t w ) = / dt'R(t,t') 



where the time enters only through the difference t — t w . 
Out of equilibrium, the generalisation takes the form 

R(t, t w ) = X(t, t w )P^-C{t, t w ) , (2) 
ot w 



(4) 



where the perturbation field is applied between time t w 
and t, that is a zero-field-cooled (ZFC) scenario. To ex- 
tract information from the susceptibility on the FDT ra- 
tio, one usually plots the susceptibility versus the corre- 
lation function and defines from the asymptotic slope of 
the curve, a number 



defining the factor X(t,t w ), the so called fluctuation- 
dissipation ratio which measures the violation of the FDT. 



(5) 
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which is often identified with the FDT ratio A^. The idea 
behind that is coming from the analyses of infinite-range 
glassy systems 22 where asymptotically, the fluctuation- 
dissipation ratio depends on time only through the corre- 
lation function. One has for those systems 

X(t,t w ) = X(C(t,t w )) . (6) 

In this context, the limiting ratio Aqo was interpreted as 
a temperature ratio between the actual inverse tempera- 
ture (3 and an effective inverse temperature seen by the 
system: [T^] /3 e // = (3X. From the dependence ©, on ob- 
tains for the susceptibility the expression 

X(t,t*)=0 f dC'X(C). (7) 

Jc(t,t w ) 

This last equation strictly holds if the violation ratio is a 
function of C(t, t w ) only and in this case, one can identify 
Aqo with A*,. It is in particular the case at equilibrium 
since X = 1 is a pure constant and it gives x = — C). 

In a recent work, |15| an exact expression was derived 
for Glauber-like dynamics which enables to calculate di- 
rectly the linear response function to an infinitesimal field. 
The advantage of this approach is obvious since it gives 
direct access to the linear response, non-linear effects are 
avoided, and not only to the susceptibility as in the ZFC 
scenario. 

In this context, we have performed a Monte Carlo 
study of the nonequilibrium evolution of the two- 
dimensional classical XY system. We have studied the evo- 
lution of the system after a quench from infinite tempera- 
ture toward the low temperature critical phase, up to the 
Kosterlitz-Thouless point. We have also considered the re- 
laxation from a completely ordered initial state. For that 
purpose, we have calculated two-point correlation func- 
tions, susceptibilities and response functions. From these 
data, we have checked the violation of the FDT and com- 
pared our numerics with theoretical predictions (spin wave 
approximation) and previous numerical works when avail- 
able. The paper is organised as follows: in the next section, 
we present the dynamics of the model and its solution for 
two-time quantities in the spin wave approximation. The 
following section deals with the numerical analyses for the 
ordered initial state. We turn then to the infinite tem- 
perature initial condition. We summarise and discuss our 
results in the last section. 



2 Two-dimensional dynamical XY model 

The two-dimensional ferromagnetic XY model is defined 
via the Hamiltonian 

H = -^2 Si- Si (8) 

<ij> 

where the sum is over nearest neighbour pairs ij on a 
square lattice and where the classical spin variables Si are 
two-dimensional vector fields of unit length. Introducing 



angular variables, one can rewrite the original Hamilto- 
nian in the form 

H = -J2 °° s (^ ~ e o)- ( 9 ) 

<ij> 

The equilibrium properties of this model are well-known 
since the pionneering work of Berezinskii. |lfi| Kosterlitz 
and Thouless and others. [T7| At a temperature Tjct, the 
system undergoes a continuous topological transition due 
to the pairing of vortex and anti- vortex excitations. Below 
the transition temperature, the system is characterised by 
a line of critical points reflecting a quasi-long range or- 
dered phase with algebraic correlation functions. The spin- 
spin correlation critical exponent rj is continuously varying 
with the temperature field. Although for the spin-spin r\ 
exponent spin-wave approximation gives an accurate an- 
alytic prediction at low temperature,^] only numerical 
estimates are known for the full temperature regime. ^Hj 
The dynamics of the model was studied extensively in 
the context of coarsening. 19 The two-time spin-spin au- 
tocorrelation function and the associated linear response 
function have been studied only recently in ref . |2()j . In the 
spin-wave approximation, valid at low temperature(T <C 
Tkt — 0.89), the nonconserved dynamics of the angular 
variable is given by the Langevin eauation|19| 

p( X , t ) = - 5 -^p-+ax,t) (io) 

where ^(x, t) is a Gaussian thermal noise with variance 
(((x, t)((x',t')) = 2TS(x - x')S(t - 1') and the free energy 
functional is given bv|19| 

F(0) = J d 2 x[V9] 2 (11) 

where p(T) is the spin-wave stiffness, related to the r](T) 
exponent by the relation 2irp(T) — T/r](T). 

Taking as initial condition a completely ordered state 
6(x,0) = 8q, using the previously defined spin-wave func- 
tional, it is possible to obtain analytical expressions for the 
two-time autocorrelation and response functions. These 
reads, at enough long times, for the autocorrelation func- 
tion C(t, t w ) = V- 1 J d 2 x(cos[9(x, t) - 0(x, t w )]):$R 

cm, 1 fa±^y' TW (12) 

where t w is the waiting time, t is the total time and A = 
t/t w is the scaling ratio. This behaviour can be explained 
in the following way. At short time difference t — t w <C t w , 
the fluctuations of small wavelength (<§; £(t w )) have equi- 
librated and we are in a quasi-equilibrium regime with a 
correlation function decaying as C(t, t w ) ~ (t — tyj)^ 11 ^^ 
where the dynamical exponent z — 2 for the 2d XY model. 
At longer times, when the scaling function significantly 
differs from 1, the aging process takes place giving rise to 
a full two-time dependence, that is a breakdown of time- 
translation invariance. The conjugate response function, 
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denned by R(t, t u 



V 



-l r ,j2 

J u x Sh(x.t 



h=0 



R[t, t w ) 



2 V (T) C(t,t w ) 

T i tin 



(13) 



with C'(t,t w ) given in equation l|12l) . It is amazing to 
notice at this point that the last equation has exactly 
the form obtained from the Fluctuation-Dissipation the- 
orem with a power law equilibrium correlation function 
C(t,t w ) ~ A(t — t w )~ v / z , where z is the dynamical ex- 
ponent. The difference from the equilibrium situation lies 
in the fact that the nonequilibrium amplitude A is ac- 
tually depending on time too. So, by differentiation, one 
has also a term arising from the derivative of the ampli- 
tude A(t w ), leading to a deviation from the FDT. From 
equations l|12fl and i|13|) . it is easy together with the defi- 
nition of the fluctuation-dissipation ratio given previously 
to obtain 



X{t,t w ) = 1 



(A-l) 2 
2(1 + A) 



(14) 



For a quench from an infinite temperature state to a 
temperature T < Tkt, no such analytical expressions are 
available. However, on the basis of scaling arguments, ^Hj 
one can postulate the general expressions 



and 



C(i, t w ) — 



R(t, ty 



1 



ft 



( m 



h 



(15) 



(16) 



where fc and Jr are the scaling function and where the 
correlation length £ has a different behaviour if the quench 
is done from infinite temperature or from a completely 
ordered initial state, namely one has: 21 



ffi 



t l/2 



Ti < Tkt 
Ti > T KT 



(17) 



The logarithmic correction in the disordered initial state 
case is due to the slowing down of the coarsening caused 
by the presence of free vortices, [21] since the approach 
toward equilibrium proceeds through the annihilation of 
vortex-antivortex pairs, a process which is slower than the 
equilibration of spin waves. 

We shall concentrate first on checking numerically 
these analytical predictions, testing at the same time the 
validity of our numerics, and then turn to the numerical 
study of the infinite-temperature initial condition. 



3 Numerics 

3.1 Numerical approach 

During the simulations, the system is initially prepared in 
two different ways: spin angles 9i chosen at random in the 



interval [0, 27r] corresponding to the infinite temperature 
initial state and constant initial angles, 9i = est. Vi, corre- 
sponding to the zero temperature initial state. For the nu- 
merical analysis, we use a standard metropolis dynamics 
where a spin chosen at random is turned at random with 
an acceptance probability given by min[l, exp(— AE/T)] 
where AE is the difference energy between the actual con- 
figuration and the former one. As stated before, in order 
to go beyond the susceptibility and to access directly the 
response itself, we use a different dynamics, Glauber-like, 
where the transition probabilities of a configuration with 
a spin Si to a new value S[ is given by 



P(Si - = 



W(S>) 



W(Si) + wis,) 



with 



m.v,): ,-xp I --s^Sj 



(18) 



(19) 



Both dynamics have the same dynamical exponents and 
one expect no significant changes for thermodynamical 
quantities. 

The two-times autocorrelation function is defined by 
C[t,t w ) - -^]T (cos [0,(0 -9i{t w )\ ) (20) 

i 

where (.) is the average over the thermal histories. In the 
metropolis simulation, we calculate the ZFC susceptibility 
via 22 



L 2 h 2 



(21) 



where h is a small bimodal random magnetic field applied 
from t w . The overline means an average over the field re- 
alizations. Practically in our simulations we use the value 
h = 0.04.E3 

The response function itself is obtained numerically 
with the help of the Glauber-like dynamics. 15 By defi- 
nition, the autoresponse to an infinitesimal magnetic field 
applied at t w is given by 



R{t, t u 



SSi(t) 



dh^tw) 

With the help of the master equation 

P({d'}, t + 1) = -> {6'})P({9}, t) 

{6} 



(22) 



(23) 



and following the lines of ref . ^H] , it is easy to arrive at 

R(t,t w ) = p(coa0i(t) [cosf^i™ + 1) -cos9f(t w + 1)]) 
+ 0{ sm6 t (t) [sm9i(t w + 1) - sin^(^ + 1)]) (24) 

where cos 9f(t) and sin 9f (t) are the components of the 
Weiss magnetisation given by 



St 



1 d 



(3 dh x 



In Z, 



(25) 



h=0 
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Zi = exp(-{3H(9i, ft)) + exp(- ( Si?(^, ft)) is the local par- 
tition function in the field. 

The thermodynamical quantities are calculated on 
square samples with periodic boundary conditions of lin- 
ear size up to L = 512 and averaged typically over 1000 
thermal histories. 



3.2 Ordered initial state 

We start with a completely ordered state and set the tem- 
perature T < TjcTt m order first to check the compati- 
bility of our numerics with the analytical predictions in 
the spin- wave approximation. In figure 1, we present the 
results obtained for the autocorrelation function in the 
asymptotic regime t — t w 3> t w 3> 1 for a temperature of 
T = 0.3 where the expression <|12|) is expected to hold. For 
different waiting times, the collapse of the data is fairly 
good and the power-law behaviour in terms of the variable 
(1 + A) 2 / (4 A) gives a very good agreement with the XY 
rj(T) exponent, as it can be seen on figure 2. 




1 10 100 

{X+lflAX 

Fig. 1. Rescaled autocorrelation function at T = 0.3 for a 
system of linear size L = 512 and for different waiting times. 
The solid line is guide for the eyes corresponding to the value 
77(0.3) ~ 0.015. 

In the asymptotic regime, the two-times response func- 
tion is expected to be given, at least at low temperature by 
the spin- wave approximation formula (|13fl . In figure 3, we 
give the numerical results obtained for a final temperature 
of T = 0.3 for different waiting times. The aging part of 
the response is very small in the accessible regime and the 
deviation from a time-invariant process is very difficult to 
be seen as the collapse of the data for different waiting 
times in figure 3 attests. Nevertheless, we have plotted in 
figure 4 directly the numerical response function together 
with the analytical prediction. The superposition of both 
curves seems to validate the expected law. Although, in 
ref. |2U| this case was considered quite extensively, it was 
done only for one temperature. Here we have extended the 
results to the whole low temperature regime. 



0.3 




0.2 0.4 0.6 0.8 1 



Fig. 2. Dependence of the 77 exponent on temperature ex- 
tracted from the two-time autocorrelation function (symbols). 
The solid line corresponds to numerical values given in ref.|18|. 
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Fig. 3. Response function at T = 0.3 for different waiting 
times. 

3.3 Infinite temperature initial state 

The infinite-temperature initial state is more canonical in 
the study of coarsening and aging effects. After the quench 
in the critical phase, the correlation length will grow in 
time with a logarithmic correction due to the interaction 
of walls with free vortices as mentioned in ref.[5Tj. That is, 
£(i) ~ (t/ln/j) 1 / 2 and leading to the conjectures (|15I16J) 
for the correlation and response functions. Berthier et al 
have checked this conjecture for the correlation length only 
for one final temperature, namely T — 0.3. In figure 5 
we show the results obtained for several quench tempera- 
tures, ranging from T = 0.1 up to T — 0.7. The collapse 
of the rescaled correlation functions (t — t w ) v ^ 2 C(t,t w ) 
as function of the variable £,(t)/£,(t w ) is very satisfactory. 
From these curves, we can extract the scaling function 
fc, see equation Ijl5|l . and find the power law behaviour 
fc{x) ~ x~ K with a temperature independent exponent 
K = 1.05(10). In ref.|2nj the value k = 1.08 was given 
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10- 




1000 



Fig. 4. Response function for different waiting times for a sys- 
tem of linear size L — 512. The quench temperature is T = 0.3. 
The solid lines correspond to the analytical expression 11311 . 



which is of course compatible with our data. However, 



times, the decay of the autocorrelation function has the 
same power law dependence as in the equilibrium situ- 
ation. Otherwise, the decay is faster with a power law 
t -1 '/ 2-1 / 2 up to logarithmic corrections. 

The parametric plot of the susceptibility times the 
temperature versus the correlation function does not 
collapse for different waiting times, showing that the 
fluctuation-dissipation ratio is not a function of C only 
but rather has a dependence on both t and t w . How- 
ever, after an initial quasi-equilibrium regime, where the 
different waiting time curves are collapsing and lead to 
the equilibrium value X(t,t w ) — 1, they have a constant 
slope, A"*,, independent of t w . This number A^ corre- 
sponds, when X(t,t w ) — X(C(t,t w )) only, to the asymp- 
totic limit Xqo = limj^^oo lim^oo X(t, t w ). In figure 6, 
where we have represented A^ versus the reduced temper- 
ature T/Tkt, we clearly see a linear behaviour, starting 
at X^ = for T = and finishing at the value X*, = 1/2 
at the Kosterlitz-Thouless point. It can be noticed that 
this continuous dependence of the FDT ratio on a model 
parameter was also observed in the spherical model, pi] 




Fig. 5. Scaling plot of the two-times autocorrelation function 
for a quench from infinite temperature towards T < Tkt- The 
three different collapsed curves are obtained for T — 0.3, T = 
0.5 and T — 0.7 from top to bottom. The different waiting 
times for each temperature are t m = 100 circles, t w = 300 
diamonds and t w = 1000 crosses. The solid line corresponds to 
l/x. 

one has to be careful in this statement since the exponent 
is very close to 1. For example, if one takes as the scal- 
ing variable x~ x = instead of x = £(t)/£(t w ), 
then the scaling limit we are interested in is <C 1 
and what is seen could be very possibly the leading ex- 
pansion terms of an analytical scaling function, that is 
~ g(0) +ax~ 1 . Moreover, numerically the extrap- 
olated value g(0) seems to be very close to zero (less than 
0.01) and it is impossible to test this value in the time 
range explored in this work. The same is true for the refer- 
ence |201 - So if g(0) is not vanishing, finally at long enough 



Fig. 6. Fluctuation-dissipation ratio versus reduced temper- 
ature T /Tkt- The line is the conjectured function. 



Finally, we present the data obtained for the linear 
response function and the fluctuation-dissipation ratio 
X(t, t w ). The simulations are done on lattices of linear size 
up to L = 100 averaged over 11000 realizations in order to 
have a good enough statistics for X. In figure 7, we show 
the result obtained for a final temperature T = 0.1. Very 
similar curves are obtained for other temperatures. The 
collapse of the data for different waiting times is very good, 
which confirms the scaling conjecture (|16fl . Although the 
number of different histories we have realized is quite huge, 
the noise on the points is still important. The range of time 
used is from t — 100 up to t — 2500, which explains the 
very short window of the x-axis in figure 7. Nevertheless, 
what it is seen after an initial short time regime is a linear 
behaviour in terms of the scaling variable £ (t) /£ [t w ) , that 
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4 Summary and outlook 



o 

0.4 - 




[(t ln 2 t w )/(t w In t)] 1/2 

Fig. 7. Rescaled response function at T = 0.1 for a linear 
system size L = 100 averaged over 11000 realizations. 



0.1 




10 20 30 

(t In t w )/(t w In t) 

Fig. 8. Fluctuation-dissipation ratio at T = 0.1 as a function 
of the scaling variable t In t w /t w In t for t w = 10 (stars), t w =30 
(circles) and t w — 100 (diamonds). 



is asymptotically 

(t - t w )^R(t, t w )^A R V2 , (26) 

with an amplitude Ar slowly varying with the temper- 
ature. In figure 8, we present the fluctuation-dissipation 
ratio X(t,t w ) as a function of t h\t w /t w hit obtained nu- 
merically at the same temperature, T = 0.1, for waiting 
times t = 10,30,100. Since in the calculation of X, we 
have to do a derivative, the results obtained are much more 
noisy and it is difficult to go on very long waiting times. 
Nevertheless, we clearly see a good collapse on a master 
curve, leading to a vanishing fluctuation-dissipation ratio 
in the asymptotic limit. The same is obtained for other 
temperatures. 



We have studied numerically the non-equilibrium relax- 
ation properties of the two-dimensional XY model initially 
prepared on two distinct ways: completely ordered and 
fully disordered state. In both cases, the two-time spin au- 
tocorrelation function and the associated linear response 
function were determined. 

In the initial ordered case, we have fully confirmed the 
analytical predictions obtained in the spin wave approxi- 
mation, strictly valid at very low temperature. Neverthe- 
less, the scaling form given in equation seems to be 
valid in a wide range of temperature below the Kosterlitz- 
Thouless transition. From it, we have extracted the equi- 
librium exponent rj(T) with a quite good accuracy as seen 
in figure 2. Using the Glauber-like dynamics defined previ- 
ously, we have obtained directly the linear response func- 
tion itself. This permits a direct comparison of our data 
with the analytical expression (|13f) . With this approach, 
we avoided difficulties inherent in the use of susceptibil- 
ities, which can be affected by short-time contributions. 
Although the aging part of the response seems to be very 
small, as attested in figure 3, the plot in figure 4 shows a 
very good agreement between the numerical data and the 
analytical prediction. 

Starting with a fully disordered state, we have ex- 
tended the conjecture checked in ref. [201 f° r one partic- 
ular temperature to the whole low temperature regime. 
For temperature ranging from T — 0.1 up to T — 0.9, we 
have numerically confirmed the forms l|15|) and i|10|) of the 
correlation and response functions with a scaling variable 
given in (|17fl . As discussed previously, we found numeri- 
cally for the asymptotic behaviour of the autocorrelation 
scaling function fc, defined in (|16f) . a behaviour which is 
compatible with a purely algebraic decay with a temper- 
ature independent exponent very close to one. The linear 
response scaling function, in the time-range studied here, 
has a linear behaviour in the scaling variable 
Those asymptotic behaviours of the correlation and re- 
sponse scaling functions are supporting, up to logarithmic 
factors, the forms given by local scale-invariance theory. 0] 
Utilising the notations of one has 

C(t,t w )~t- a F c (t/t w ) (27) 
R(t,t w ) -t-^FaWtv) , (28) 

where the scaling functions Fc and Fr have the asymp- 
totic forms 

F c>R (u) ~ A c ,r vT Xc ' R l z u » 1 . (29) 

From our data, we obtain a = rj(T)/2 and Xq = Xr = 
r)(T) + 1 confirming the general scenario depicted in ref.0J 
8 ,25 . Finally, we extracted X^, from the parametric sus- 
ceptibility/correlation plot in the long-time limit. The re- 
sult we obtained is well fitted by the linear behaviour 
Xg; = (1/2)T/Tkt- The direct use of the response func- 
tion gave a different answer, as seen in figure 8. Al- 
though the fluctuation-dissipation ratio X(t, t w ) is a func- 
tion of both t and t w , this dependence seems to enter only 
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through the scaling ratio t lnt w /t w \nt. From the numer- 24. A. Picone and M. Henkel, J. Phys. A 35, 5575 (2002) 
ical results we obtained on the correlation and response 25. M. Henkel, M. Paessens and M. Pleimling, Europhys. Lett, 
functions, it is clear that in the asymptotic regime the flue- 62, 664 (2003) 
tuation dissipation ratio X x = limt M _»oo lim^oo X(t, t w ) 
is vanishing. A result which is different from what is ob- 
tained from the parametric susceptibility/correlation plot . 
This vanishing is due to the breakdown of scaling in- 
duced by the presence of the logarithmic factors in the 
scaling functions. In practice, one has to take care when 
discussing those plots, especially when no master curve is 
ever reached for different waiting times. 
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